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Abstract: A Markovian Monte Carlo algorithm for multi-parton production in the high-energy 
limit is proposed and the matching with unintegrated parton densities is discussed. 



1 Introduction 



Hard scattering at hadron colliders is usually de- 
scribed within the framework of collinear factorisation 
[TJ[2]. The full scattering amplitude is factorised into a 
hard perturbative parton scattering matrix element and 
process-independent universal parton distribution func- 
tions (PDFs), which depend on the flavour of the ex- 
tracted parton, its energy or light-cone momentum frac- 
tion x w.r.t. the initial hadron, and on the factorisation 
scale (j, p. The choice of this scale is, to some extent, ar- 
bitrary. By the inclusion of higher-order corrections the 
dependence of the cross section on this scale is diminished. 
At present, PDFs cannot be obtained from first principles 
due to their essentially non-perturbative origin, but they 
can be extracted from data, for example through global 
fits [SI HI IS]- On the other hand, the evolution of these 
collinear PDFs with changing factorisation scale can be 
determined perturbatively. In the collinear factorisation 
scheme, all initial state partons are on-shell and have zero 
transverse momenta k_L = 0. 

An alternative approach is the framework of kj_- or high- 
energy factorisation. There, unintegrated PDFs (UPDFs) 
are convoluted with off-shell matrix elements. The PDFs 
are unintegrated in terms of the initial partons' kj_. Ini- 
tially, kj^-factorisation has been formulated for heavy 
quark production J6][7j[8]. The approach has been further 
investigated in other channels, see for instance [51 HD1ITT] . 

The k^-factorisation has apparent advantages over con- 
ventional collinear factorisation: First, in the high-energy 
limit, i.e. for t <C s with s being large, the QCD cross 
section for jet production is dominated by gluon exchange 
diagrams, which diverge in this limit. This divergence is 
alleviated or even removed by realising that the 1/t diver- 
gences in the matrix element can be identified with diver- 
gences of the form l/kj_ and thus using a suitable form 
of unintegrated PDFs, vanishing fast enough for k^ — * 0. 
Second, employing UPDFs means including the leading 



logarithmic contribution of higher order corrections to a 
given process, since the effect of additional QCD radiation 
is encoded in them [T21H5] . 

Taking the high-energy limit in a given process is equiva- 
lent to the BFKL limit [T4"lll5j . which builds on i-channel 
dominance of scattering cross sections and the reggeisation 
of i-channel gluons [15]. In the past, there have been vari- 
ous approaches, aiming at a solution of the BFKL dynam- 
ics with Monte Carlo methods and thus producing exclu- 
sive final states. An approximation, aiming at a correct de- 
scription of essential features of the BFKL equation and a 
correct extrapolation to the DGLAP regime, has been pro- 
posed in the "Linked Dipole Chain Model" [171115] . This 
model has been implemented in [19 . The scope of this ap- 
proach is closely related to the CCFM equation [2T?ll2ip22| . 
Event generators based on this evolution equation have 
been presented in 23,24,25,26]. An iterative solution of 
the pure BFKL equation has been proposed in [57] , itera- 
tive Monte Carlo solutions in [2"8 ] |29 ] . Later on, this pre- 
scription has been extended to next-to-leading logarithmic 
accuracy (30p31 l l3"2" l l3"3"]. 

In this paper, a different implementation of k^-factorisa- 
tion for the case of multijet production is discussed. Em- 
phasis is put on finding a gauge invariant form of the cor- 
responding expressions and on identifying their matching 
to unintegrated PDFs derived from conventional collinear 
ones [SHISSIISS]- It turns out that this in fact can be 
achieved by working in the high-energy limit, using as ba- 
sic building blocks splitting functions in the limit z — ► 00 
in conjunction with a proper reggeisation of all i-channel 
propagators. Since four- momentum conservation can ex- 
plicitly be imposed in a Monte Carlo solution, this ap- 
proach clearly includes effects beyond the naive leading 
order BFKL limit @ Furthermore, identifying the proba- 



1 In addition to the pure gluonic ladders of the high-energy limit, 
here also vertices for quark production are included. 

2 As was discussed for example in 37,38,39 32 , the implementa- 
tion of four-momentum conservation and running a a effects strongly 
modifies naive LO BFKL predictions, which were shown to poorly 
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bilistic interpretation of each emission in the high-energy 
limit, the Monte Carlo solution has for the first time been 
implemented as a Markovian approach, similar to conven- 
tional parton shower event generators. This enables gener- 
ation of an a priori arbitrary number of emissions, which is 
important at high energies, where corrections due to large 
final state multiplicities are sizable. 

The paper is organised as follows: In Sec.[2]the procedure 
of [34,35,36,40] (KMRW) to generate doubly unintegrated 
PDFs (DUPDFs) and the corresponding angular ordering 
constraints are reviewed. In Sec. [3J it is then shown that 
the leading ln(l/x) terms are correctly taken into account. 
Section H] contains the description of the Markovian MC 
procedure to generate event topologies with an a priori 
undetermined number of final state partons. In Sec. [5] first 
results are presented and Sec. [7] contains our conclusions. 



2 Unintegrated parton densities and 
the KMRW procedure 

In this section, the KMRW procedure of constructing un- 
integrated PDFs from conventional DGLAP PDFs PPS] 
36, 40] is reviewed. The discussion and notation closely 
follows [361140]. 

In collinear factorisation, where the transverse momenta 
kj^ of incoming partons are taken to be zero, the par- 
ton densities obey the DGLAP evolution equation [HJ 
H21|43l[44] , which determines the /^-dependence at fixed 
light-cone momenta. This evolution equation resums lead- 
ing logarithmic parts of higher perturbative orders. In a 
Monte Carlo formulation, real emission corrections can be 
implemented as a Markov chain of 1 — ► 2 parton split- 
tings [45,46,47 ]. However, a study of QCD beyond dou- 
ble leading logarithmic order reveals that quantum coher- 
ence effects suppress parton emissions in regions of phase 
space, where the emission angle of the emitted parton 
is larger than the opening angle of the emitting colour 
dipole [48,46 . To exemplify this, consider a parton evolu- 
tion chain in the initial state of a DIS event, as depicted in 
Fig. [TJ If angular ordering is fulfilled, the momenta k% of 
the radiated partons will be distributed such that their an- 
gle 8i with respect to the beam direction increases from the 
incoming proton towards the hard scattering. To investi- 
gate the implications of this constraint, it is convenient to 
start with a Sudakov decomposition of the momenta [55] , 



Pi 



x,P + ftg' - k, 



k i = a i P + /3 i q' + k i x, (1) 



where P is the proton momentum, q is the photon momen- 
tum and q' = q + xbP, with xb being the Bj0rken x. In 
the high-energy limit, the proton mass can be neglected, 
m 2 <C Q 2 = —q 2 - Hence q' 2 = and in the Breit frame 




Fig. 1 Multiple gluon emission in deep inelastic 
lepton-nucleon scattering. The hard scat- 
tering process is characterised by the scale 
fi. Usually this scale is also employed as the 
factorisation scale. 



the momenta read 
1 



P 



(0,0, Q) , 



2x B 

q' = l(Q,0,-Q) and 



ki± = (0,ki_L,0) . 

All emitted partons are on-shell, which allows to relate 
their Sudakov parameters through 



ft = (ft_! - ft) 



I - Zi Xijx B 



where z% = Xi/xi—i. Imposing angular ordering for the 
emissions results in ordering of the corresponding rapidi- 
ties yi, since 

Vi = ^ m 6 = - In tan ^ , 

where = kf /k^~ — oh/xb ft and 9i is the angle of ki 
with respect to the beam axis. According to Eq. |T]) 



£ ^ i ( ^ %i Q ^\ x^ ( Q 



Xb V %i kj 



(2) 



where the rescaled transverse momentum kj = kj Zi) 
has been introduced. Hence angular ordering requirements 
yield the constraints 



Ziki < k i+1 and z n k n < p 



(3) 



describe data. 



Herep = x n+ iQ\^./xB is the maximal rescaled transverse 
momentum which is fixed by the hard process through 
S = (1 + (3 n +i)/(x n +i/xB - 1)- Typically, in an angular 
ordered evolution of the parton distributions, p plays the 
role of the factorisation scale \ip [2Tni2"TU2"2l5"U] . The above 
ordering procedure can be generalised to hadron-hadron 
collisions. In this case, both incoming particles have a 
partonic substructure. In general, this leads to two sepa- 
rate factorisation scales, yt^) and fj,p\ for the two parton 
densities, respectively. 

In [34, 35, 36,110] it has been shown that doubly uninte- 
grated PDFs (DUPDFs) may be inferred from conven- 
tional DGLAP PDFs. In the following, DUPDFs will be 
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denoted by f a (x, z, k\,H F ), while their conventional DG- 
LAP counterpart will be denoted by f a (x,Hp)- The DU- 
PDFs must satisfy the normalisation condition 

f 1 f dk 2 

J dz J -p- f a (x,z,k ± ,H F ) = xf a [x,fi F ). (4) 

Employing the Sudakov form factoJl 
A„(ki,^) (5) 



exp 



dkl a s (k'l) 1 
2 k /2 i 2?r 2 ^r 1 



d(P ab (0 



with P a b{C) denoting regularised DGLAP splitting func- 
tions for the splitting a — ► 6, a singly unintegrated parton 
distribution / a (x, k^, /Li^) is obtained through 



f a (x,k\, n 2 F ) = 



d 



xf a {x,k\)h a {k\,H F ) . (6) 



dhxk\ 

In the region k\ > Hf this UPDF is set to zero. This 
procedure leaves some minimum k^-scale to be defined, 
below which DGLAP parton evolution is not valid. In the 
following, this scale will be denoted by /it 2 ,. Relation © 
then holds true only above Ho, which yields the constraint 



^ dkj 
k 2 



fa(x,k\,iJp) = xf a (x,nl)A a (ij,l,Hp) 



on the singly unintegrated PDF. Whenever UPDFs, sat- 
isfying this normalisation condition, are applied in kj_- 
factorisation, physical observables must be insensitive to 
details of the infrared behaviour of f a (x, k 2 , up), i.e. be- 
low //qO Therefore, a choice can be made, for example [ID] 



fa(x, z, kj_, Hp) 



Mo 



2 



k 2 



A a {x, z, up) H {x, z, /ip) 

Mo 



where 



A a (x, z, Hp) — — fa (x, z, Ho, Hf) 
2x 

+ Y~^. f a ( X > Mo) A aOo; Mf) 

B a (x, z, hf) = 2 f a (x, z, hI, Hf) 



2.T 



1 



fa(x,Ho) ^o(MO'Mf) 



This choice implies that the UPDF vanishes with k^ for 
k± — » 0, as required by gauge invariance [51] . 



3 The factor of 1/2 in the sum over the parton species avoids 
double-counting s- and t-channel partons. 

4 It turns out that there is no need for an explicit form of the 
DUPDFs below /Xg, since the t-channel parton chains contain a nat- 
ural cutoff in k^, cf. |28] . by imposing phase space cuts given by 
physical observables like minijets. 



Instead of the regularised splitting functions P a b(z), un- 
regularised splitting functions P a b(z) may safely be used 
here. This is because the splitting kernels are implicitly 
regularised by imposing the rapidity ordering constraint 
Eq. ([3]). Inserting corresponding O-functions in z results 
in the singly unintegrated quark and gluon distributions 
f q (x,k 2 ± , Hp) and f g (x,k 2 L , Hp), respectively pOTTO] , The 
term singly unintegrated indicates that these PDFs de- 
pend on one additional variable w.r.t. the collinear ones. 
It is straightforward, however, to introduce an additional 
z-dependence by simply dropping the z-integration in 
Eq. ©. Such defining the DUPDF 



f a {x, z, k 2 ± , h 2 f) 



A a (ki, Hp) 



(7) 



(1 - 5 ab ) + S ab Q 



Hp 



Hf + kj 



the desired relation, Eq. |4]), is satisfied for both parton 
species. To guarantee the consistency of the approach, 
the conventional DGLAP PDF employed to obtain the 
DUPDFs should be determined using the leading order 
unregularised splitting kernels employed in Eq. JJJ). Fur- 
thermore, a consistent treatment of the running coupling 
a s should be imposed. 



3 DUPDFs as impact factors for LL 
BFKL evolution 

In this section we argue that the DUPDFs defined above 
may be employed as impact factors in the calculation of 
multi-gluon cross sections in the high-energy limit. The 
argument works at leading logarithmic (LL) accuracy. The 
starting point is the integrated LL gluon branching prob- 



ability r 



LL) 



= -logA^ 



which determines the be- 



haviour of the DGLAP evolution of the gluon density^ 

r("V,M 2 ) = r(f V,m 2 ) + £i^ l V,m 2 ) , 

where 

rlt L V,M 2 ) = [ ln " 2 dlnki f^dz^PM 



with P ab (z) again denoting the unregularised DGLAP 
splitting kernels and the integration boundaries deter- 
mined by angular ordering, cf. the 0-function in Eq. ([7|. 
To simplify the discussion we firstly focus on Tg^ L ^ only. 
The corresponding part of the Sudakov form factor reads 

A(f)( A1 2 ,M 2 ) = ex P {-r(f)( A1 2 ,M 2 )} . 



5 The factor 1/2 contained in Eq. (0 must be cancelled here in 
order to restore the t/n-symmetry of the splitting process. 
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Replacing the splitting variable z of the emitter parton by 
the rapidity y of the emission, which, according to Eq. ^ 
is given by 



y = o m £ = ln v i 

2 \xb k 



results in 



In- 



In /I 2 



dy 



dlnki 

ln/i 2 ^J/(2min) 

2C A (l-z{l~z)f a s 

P gg { Z ) 2vT 



(8) 



dlnki / d y a s (l-z(l-z)) 2 , 

ln/1 2 ^y(Zmax) 



where a s = UsCaI'R- The term z(l — z) in the numera- 
tor corresponds to helicity non- conserving configurations 
in the 1 — > 2 parton splittings and thus in the impact 
factor [52] ■ These configurations are absent in the high- 
energy limit, which simplifies the integrand of Eq. (|8j). 
such that the part of the integrated LL gluon branching 
probability induced by g — > gg splittings reads 



rlnft 2 /-J/(zmin) 

Tf g L Hp 2 , ft 2 ) = / dlnki / dya s 

J In /j- 2 J y(2 raax ) 



(9) 



Keeping in mind that a s depends on transverse degrees of 
freedom only, performing the y-integration results in 



l - 2 

• In ft 



dlnk^ a 



In fi 2 
X ihl 



ft xQ 
kj_ x B \tx 



-ln 



k^ xQ 
ft x B ^i 



2 Jo dln2 ki^ 



The order of integration in Eq. (J9j> may be changed, 

WA-J./V/ ^ (10) 

xe(lnA 2 /ki+y-y') , 

where y = Iux/xb + hiQ/ft and y ~ y — \nft 2 /^ 2 . 
If the running coupling is treated identically, this result 
agrees with the reggeisation factor of the i-channel gluon 
propagator found by rewriting Eq. (7) of [35]. Up to a 



minor transformation, this equation read^l 

/" {Vab, PaX, Pbx) 
. n 

n 



(ii) 



Hk 2 

a s dy 



q 2 

,o — exp<j -a s ln-^Ay. ; 

kf 2tt I ^ 



x exp 



q 2 1 1 
-a s ln -^-Ay -5 (jp bl _ + q n ±) 
Mo ) z 



where a s — a s G'A/ir and qi — p a — J2]=i kj- The expo- 
nential term in the square brackets is readily identified as 



A(y,y) = exp{-f( iL )(y,y)} 



(12) 



where 



tf L \y,y) = [" dy 1 f 

Jv JO 



d ln a s , 



which is the desired result. It has been pointed out e.g. 
m [T3] that the comparison with NLO BFKL calculations 
suggests the choice a s = a s (k 2 L ), similar to the DGLAP 
case. Employing 



f3 logk 2 JA 2 ' 



where /3q 



11 - 2/ZNf 



-In 



we then end up with the result presented in 



n LL) (y,y) = (y-y) % log {^§{ 



In our numerical analyses, A is chosen consistent with the 
input PDF. Equation ITT]) can be used to construct the 
full LL BFKL kernel / through 



f(Vab, PaX, Pbx) = ^2 P {Vab, PaX, Pbx) 



(13) 



n=0 



Since rapidity ordering is trivially satisfied in the BFKL 
evolution, the explicit ordering requirement incorporated 
in the 0-function of Eq. (|I0|) may be dropped whenever 
A(y, y) is employed. 



is given by 



Following the same reasoning, T g ^ L ^ 

•In /i 2 ry(Zmn) 

d ln k^ / dy a s 

lnp 2 ^J/(Zmax) 

xg^z(l-z)(z 2 + (l-z) 2 ) • 

In principle, this term vanishes in the high-energy limit 
due to the prefactor z(I — z), thus allowing to identify 
A(y,y) with A g LL \fi 2 , ft 2 ). However, it may be used to 
model quark production along the BFKL ladder, as will 
be discussed in Sec. [5] 

6 Note that the particle indices a and b are interchanged with 
respect to Schmidt's original formulation. 
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Similar considerations may be applied to the integrated 
quark branching probability. Starting from the expression 

r(")(^,A 2 ) = r(f)( M 2 ,M 2 ) + r(f)(^,A 2 ) 

and again replacing the splitting variable z by the rapidity 
y results in 

■In /i 2 /•Iz(Zmin) 

dlnki / dy 

x5 s ^ (l-z)(l + (l-z) 2 ) . 

By identifying z — —t/s, all factors 1 — z become unity in 
the high-energy limit. Thus, 

Simultaneously, due to the denominator part (1 — z) in 
P qq (z) quark production in the i-channel is suppressed, 
hence allowing to identify 

However, T q ^ L \[i 2 , jl 2 ) may be employed to model gluon 
emission from t-channel quark lines, as will be described 
in Sec. [5] 

The above considerations show that to leading logarithmic 
accuracy the DUPDFs, Eq. ([7]), resemble all features of the 
BFKL evolution. Therefore, they can safely be employed 
as impact factors for the calculation of cross sections in 
the high-energy limit. 



4 Markovian Monte Carlo solution to 

the ln(l/x)-evolution 

The Markovian approach to the calculation of cross sec- 
tions and differential distributions in the high-energy limit 
will be presented in this section. The advantage of the 
algorithm is that the number of emissions stays a pri- 
ori undetermined, similar to the case of conventional 
parton showers employed to solve log(Q 2 //x 2 )-evolution 
45,46,47 . The factorisation of the radiation pattern 
into individual emissions, which depend on each other 
merely through the correct ordering, allows to model fur- 
ther physics effects involving the produced outgoing par- 
tons, like for instance adding final state radiation. 

The basis of the formalism is encoded in Eq. (7) in [28] and 
Eq. (|12p . These equations translate into the probability 
for having an additional emission from the BFKL kernel 
being approximately distributed according to the function 



n (l,Tf L \ yi ,y n )) 

= Tf L \y^y n ) exp{-r(")(^,y n )} 



(14) 



Here, yi is the rapidity of the previous and y n is the ra- 
pidity of the final emission. Such distributions may be 
generated employing the veto algorithm, described for ex- 
ample in [47]. It allows to simultaneously select the ra- 
pidity and transverse momentum of the new emission In 
the following, the superscripts will be dropped. 

To determine the corresponding z-k^-factorisation formu- 
la, the simplest case, a gluon ladder with no emission, is 
investigated. This corresponds to a "2 — > process" in the 
z-kj_-factorisation approach. When working in collinear 
factorisation rather than with the DUPDF prescription 
of [JO], it is a 2 — * 2 process. The corresponding phase 
space element can thus be determined by factorising the 
collinear matrix element and its phase space integral. The 
starting point is 



- = E / d ^ (1) / ^ 



:(2) 



d 4 fci f d 4 fc 2 



x 6 (k 2 ) 8 (k 2 ) (2n) A 5^(P - fci - k 2 ) 

X JaW( x '.y )JaW( x ,Q ) 2^(l)g(2)g 2 



(15) 



where the factor 1/2 is due to the identity of the final 
state particles, Q 2 denotes the factorisation scale, P 2 = s, 
s = £V-'£@)S, £ = x/z, and the superscripts W and ^ 
refer to the left and right beam, respectively. The matrix 
element reads 



\M, 



mil 



(47ra s ) 



r 2 



tu us st 



(16) 



Employing Z\ = z 2 = z, t — —zs and u — —(1 — z)s 
transforms this into 

\M gg \ 2 = (4na s ) 2 ±[P gg (z)] 2 {l + 0(z(l-z))} 

where terms proportional to z(l — z) in the numerator 
vanish in the high-energy limit and are not explicitly dis- 
played. 

The phase space element of the general case of a gluon lad- 
der with an arbitrary number of gluons emitted between 
the two outermost jets can be derived by combining their 
momenta into one final state momentum K . Ignoring the 
substructure of K , the differential two-particle initial and 
final state phase space element for the remaining degrees 
of freedom reads 



d$ 2 = d£ (1) d£ 



. (2) d 4 fcr d 4 fc 2 . 2 . 
(2tt)3 (2tt)3 6 K l) 



6 (k 2 2 ) 



x (27r) 4 <5 (4) (P-K-h-h) , 

with P again the total four momentum of the process. 
Employing the four-dimensional (5-function and the rela- 
tions d^Wd^ 2 ) = dyds/S and dp z — d (y/s± sinh y) = 



7 In fact applying a veto is not necessary here, as long as quark 
production is neglected in the approach. 
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sj_ cosh ydy = Edy results in 
2tt dyi dkf ^d^i 



d<P 2 = -dsdy 4(2?r)3 



Furthermore, the definition P = P — fo allows to rewrite 



dy d 1 , P+ + m 2 ±e+y 2 
- In ■ 



dy2 dy 2 2 P~ + m 2 xe~ y2 
1 / m2j_e +J ' 2 m 2 j_e~ V2 



2 V ^ 



P- 



Pfca 



Using P = y/s (coshy, 0, sinh y) gives 
ds 5 (s + K 2 - 2P(K + kx) + 2Kk x ) 



s-P(K + h) Pk 2 
such that 



d$ 2 



1 



45 (2tt) 



dy 2 dyi dkj A 



Finally, when fixing the factorisation scale in Eq. (|15p and 
the renormalisation scale in Eq. (| 16(1 to be the transverse 
momentum in the process and adding a Regge suppres- 
sion factor for the t-channel gluon, the z-k^-factorisation 
formula reads 



25 



dyi / dk 1± / d</>i / dy 2 



x f g (xW,z,^i)f g (xW,zX ± ,tii) 



(2) 



2^(1)2^(2)25 A g (y u y 2 ) 



Here, f g is defined such that only gluon splittings are 
contained in the sum over parton species of Eq. ([7]) and 
angular ordering is implemented by the DUPDFs, while 
A g (yi,y 2 ) is given by Eq. (TT2"|) . The superscripts W and 
( 2 ) refer to the left and right beam, respectively. Since the 
emitted gluons are distinguishable due to rapidity order- 
ing, the symmetry factor 1/2 appearing in Eq. (|16p must 
be dropped. The factorisation scale hf of each DUPDF 
introduced in Eq. ([7]) is unambiguously determined by the 
rescaled transverse momentum k_L of the emissions. 

Equation (fTTjl describes a gluon ladder with no rung, 
but it can be easily extended to final states with an ar- 
bitrary number of gluons. In contrast to the previous 
case, the momentum fractions z^' and z^ are then gen- 
erally different from each other. Hence we define the 
rescaled transverse momenta kl 

K n-1± 



(1) 



, 2± - k 2± /(l - Z 

k„_ 1± /(l - zW). Employing Eq. (7) of 



and 
, the 



cross section for the 2 — * n gluon scattering reads 

a = Jsj Ayi J dk?± / d01 / dyn (18) 

X f CtW Z W k 2 k (1)2N l f (t^ Z^> k 2 k (2)2 1 
X J g [£ ,Z ,K 1± ,K 2J _ )Jg\-L ,Z , n_L' rt— IX/ 

57 £ 



2 ^(i) 2 ^(2)2 5 A g ( yi> y 2 ) [f = \J 2vr 
x [d Vi [ S k 2sMt) Clw A,( WlW _ 1 ) 



where 



C a 



C A 



The corresponding Monte Carlo event generation algo- 
rithm can be described as follows: 

1. Determine the kinematics of the first emission and the 
rapidity of the last emission according to the modified 
z-fc^-factorisation formula, Eq. (fT8|) . 

2. As long as phase space allows, choose a new rapidity 
2/i according to Eq. (I14[) and a new transverse momen- 
tum kjj_. The corresponding cuts on the individual 
emissions have already been discussed in |28j . In the 
notation employed ibidem, they are given by 



k li > Mo and 



3. Fix the transverse momentum of the last emission 
through overall momentum conservation. 



(17) 5 Model for quark production 



So far, it has been shown that Eq. (fT8|) yields the correct 
LL gluon evolution in the high-energy limit. In this limit 
quark production is strongly suppressed due to the spin 
structure entering the corresponding vertices. However, 
energies and rapidity intervals at real colliders are finite 
and quarks do appear as final state partons. Since, for 
instance, heavy quark production is of large phenomeno- 
logical interest, it needs to be described. In our approach 
we aim at not spoiling the high-energy gluon evolution. 
Therefore we choose to model quark production within the 
BFKL ladder structure by simply adding a g*q* — > q effec- 
tive vertex, which vanishes in the high-energy limit, but 
keeping the finite, non-leading terms. Additionally, quarks 
can be produced by employing qg* — > q and qq* — > g 
impact factors contained within the DUPDFs. These 
quarks may further radiate gluons, which is modelled by a 
q*q* —>■ g vertex. Figure [5] shows a possible configuration 
of quark production. 

Following Sec. El the remaining vertices are then read- 
ily determined. At leading logarithmic accuracy they are 



G 







Pa 




k n 



o 



qi-i 



pb 



Fig. 2 Multi-Regge amplitude including the emis- 
sion of a quark pair with the particle indices 
i and i + 1. The shaded blobs represent the 
vertices proposed in Eq. (1191) , 



given by the corresponding DGLAP splitting functions in 
the high-energy limit, 



Cqg — Cf , 



2 



(19) 



Then, the general case of a parton cascade in the high- 
energy limit reads 

a = |^E/ d yi/ dk ^/ d ^/ d ^ ( 2 °) 



^Coi-.ai (qi-i,ki) A ai (yi,yi_i) 



where now both quarks and gluons are contained in the 
sums over parton species. 

If heavy quarks are included in the simulation, their 
masses are taken care of in the Reggeisation factor and the 
phase space integration. Following the discussion in 53 1, 
the branching probability Tq 1 (y, y) for heavy quarks of 
mass m is modified by 



C' 



k 2 ± + m 



2 Cgg 



(21) 



Accordingly all external momenta are constructed employ- 
ing the correct on-shell masses of the corresponding par- 
ticles. 



6 Results 

In this section, results obtained with the Monte Carlo al- 
gorithm described above will be presented. All of them 
have been obtained with an implementation into the 
MC event generator Sherpa 54j|3 To eliminate possi- 
ble dependencies on the phase space integration, we have 



rp 0.6 

> 
a 
C5 



0.4 



0.2 



running « s ( x 0.1 ) 
fixed a s 

fixed a s , analytic 



k 1j= 50 GeV 
a y = y -y„ = 4 
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Fig. 3 
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k„I GeV ] 

Transverse momentum spectra / (k„x) for fixed 
and running coupling solution of Eq. (I18|l at fixed 
kij_ and Ay. Note that the result for running cou- 
pling has been rescaled by a factor of 0.1. 



cross-checked our calculations with a different integration 
method. This method uses an iterative approach to gener- 
ate event topologies for a fixed number of final state parti- 
cles, as explained in the appendix. We found no deviations 
from our results generated in the Markovian approach. 

Firstly, we focus on purely gluonic processes, reflecting 
the behaviour of the LO BFKL equation. This essen- 
tially translates into invoking Eq. (|18p for event gener- 
ation. In Fig. [3] the azimuthally averaged k„j_ spectrum 
/ (k„jj = (/ (k n ±)) < p is shown, where we have fixed ki^ = 
50 GeV and Ay = 4, and where the DUPDFs have been 
set to 1. Therefore, this plot investigates the behaviour 
of the BFKL kernel, Eq. (flU)) . only. As collider setup, 
the LHC with a cm. energy of 14 TeV has been chosen. 
In the fixed coupling solution a s has been evaluated at 
scale k1 ± - The figure shows the effect of going from a 
fixed coupling and unconstrained kinematics to a running 
coupling with kinematical constraints, which considerably 
widens the distribution. Also, since a s is typically eval- 
uated at smaller scales, / is significantly enhanced. The 
large influence of kinematical constraints and running cou- 
pling on the BFKL dynamics has already been noted, e.g. 
in [321I39]. 

As a next step, jet-production is investigated, compar- 
ing the results of the new algorithm to those obtained 
in collinear factorisation with on-shell matrix elements, 
which in the following will be denoted by DGLAP. The 
DGLAP results have been subject to the following correc- 
tions and constraints: 



• ordering of final state momenta in rapidity, 

• setting /i^ 2 = kf ± and pip 2 = k^ ± , 

• evaluating the coupling weight as Yii a s (k|i) 



'This code is available from the authors upon request. 
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Fig. 4 



3 3.5 

log ( kj/ GeV ) 

Comparison of log (kx)-distributions between 
BFKL and reweighted DGLAP matrix elements. 
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Fig. 6 Comparison of kx-distributions between BFKL re- 
sults with and without the inclusion of quarks in 
the simulation. 



However, without any i-channel reggeisation factor in the 
DGLAP matrix elements there are still large differences. 
Applying a i-channel reggeisation weight to the DGLAP 
calculation results in much smaller discrepancies. The cor- 
responding comparison for the log (kj_)- and Ay-spectra is 
shown in Figs. Q] and [5] Due to the formal equivalence of 
Eqs. (fPof and (fT7|) at leading logarithmic accuracy, agree- 




0123456789 10 



Ay 

Fig. 5 Comparison between BFKL and reweighted 
DGLAP matrix elements for the Ay-distributions. 



ment is to be expected and can be interpreted as another 
indication for the validity of the approach. Sizable devi- 
ations occur for k± > 5 GeV, which is due to the fact 
that the BFKL approach is bound to describe large en- 
ergy partons only incompletely. In order to verify this, 
we have reweighted the BFKL matrix elements with the 
exact matrix element obtained in collincar factorisation. 
The corresponding correction weight for & 2 —> n gluonic 




0.5 1 1.5 2 2.5 3 

log(k bl /GeV) 



Fig. 7 Comparison of log (kbx)-spectra, calculated ei- 
ther using the matrix element given in [H] con- 
voluted with DUPDFs or employing Eqs. ([20]) 
and l2l]). 
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k'f [ GeV ] 



Fig. 8 Comparison of jet-k^-spectra with CDF data. De- 
tails of the analysis can be found in [55] . Dashed 
lines show contributions from subsamples of 2- to 
4-particle final states. 




kf [ GeV ] 



Fig. 9 Relative differences in jet-kx-spectra compared be- 
tween the Monte Carlo results and CDF data in 
Fig. El 



process reads 

_ 8n!M gg _ >TOg (l, 

W_ (4 7 r as ) 2 F g3 (zW)F 33 (z(2))nr = " 2 1 16^ 2 a s /k^ ' 

where the factor n\ occurs due to the rapidity ordering 
in the BFKL approach and cancels the symmetrisation 
of the full DGLAP matrix element A/ gg _> ng . Performing 
this reweighting yields exact agreement between the two 
approaches. 

In a next step, all possible parton splittings in the 
DUPDFs as well as in the BFKL kernel have been en- 
abled, i.e. Eq. PU|) has been employed. As can be seen in 
Fig.® this results in a significant change of the k^-spectra 
of the partons in the high-k^ region, which is mainly due 
to the fact that quarks from the PDFs tend to have larger 
energies than the gluons. To examine the additional ef- 
fect of heavy quark masses, we have compared our results 
to those obtained in high-energy factorisation along the 
lines of [6]. For this comparison we have used the full 
off-shell matrix element convoluted with DUPDFs. Fig- 
ure [7] shows the log (kjj-spectra of the heavy quarks in 
6&-production. The coupling weight in the matrix element 
of the high-energy factorisation approach has been set to 
a s (k^j_) a s (k| ± ) in order to match the coupling weight 
in our approach. We obtain reasonable agreement with 
our calculation for k^ > 2m;,, where mass effects beyond 
Eq. pTj) are expected to have less impact on the results. 

Finally we have compared our results to recent experimen- 
tal data. Firstly we show a comparison to data obtained 
by the CDF collaboration [55 . The corresponding pre- 



diction of jet-kj^-spectra from our MC implementation is 
shown in Fig. [HI It fits the data considerably well, both 
in their shape and their normalisation. Note that no K- 
factor has been employed in the calculations. Although we 
observe a tilt of the distribution, which potentially arises 
from missing s-channel contributions to quark production, 
this is a quite remarkable result considering the fact that 
we employ a modified LO BFKL kernel for event genera- 
tion. As can be seen in Fig. [9l deviations are up to «50%, 
which is well within the expected leading logarithmic ac- 
curacy. 

Secondly we compare to the decorrelation observable in- 
vestigated in Ref. 56J . As can be seen in Fig. [10] our ap- 
proach does not completely describe the data. However, 
the deviations are of similar size as in Ref. . We stress 
that the data have not been corrected to the parton level 
and therefore correlated and systematic errors might have 
an impact. 

7 Conclusions 

In this publication we have presented a new Monte 
Carlo algorithm for the description of particle produc- 
tion through the BFKL evolution equation. This has 
been achieved in a Markovian approach, iterating inde- 
pendent emissions in order to obtain the full BFKL radi- 
ation picture. It has been discussed how doubly uninte- 
grated PDFs, obtained from conventional PDFs through 
the KMRW procedure can be employed as impact factors, 
retaining essential features of small- a; physics encoded in 
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Fig. 10 Comparison of the jet decorrelation observable 
presented in [56] with DO data. The full error 
bars include both statistical and systematic er- 
rors, whereas statistical errors are independently 
highlighted by the smaller error bars. 



the BFKL equation. In our opinion, this constitutes an 
important step towards a more unified event description, 
which allows to employ conventional PDFs deduced from 
global fits, rather than specialised parton distributions. 

The implementation of this algorithm within the frame- 
work of a multi-purpose event generator has begun, and 
first results have been discussed. They indicate that the 
proposed algorithm correctly reproduces the BFKL fea- 
tures visible in analytical calculations as well as in other 
MC approaches. The results also show the important 
effect of a running of the coupling and of kinematical 
constraints, which go beyond the LO BFKL approach. 
The realisation of the Markovian algorithm is comparably 
straightforward. Using DUPDFs obtained from collinear 
PDFs allows to compare our results for jet production to 
those obtained in the collinear factorisation approach. We 
found that we can obtain good agreement between both 
approaches, even for multi-parton production, when ef- 
fects that are not present in both approaches, such as t- 
channel reggeisation and rapidity ordering, are taken into 
account. In the same framework a model for quark pro- 
duction, which is beyond the LL approximation, has been 
included and its effect on jet production has been stud- 
ied. Finally, we found that the new approach is capable 
to describe the production of high-kj^ jets at the Tevatron. 

This work is a first step towards a unified description of 
particle production in the regime of high and low trans- 
verse momenta, i.e. of jet- and minijet-production. The 
formalism presented here can be extended to the simula- 
tion of multiple parton interactions, which constitute an 
important part of the underlying event. Also diffractive 
processes and quarkonia production may be included in 
the description. 



Acknowledgements 

We like to thank A. D. Martin and M. G. Ryskin for fruit- 
ful discussions. We are especially grateful to J. R. An- 
dersen for discussions concerning the treatment of run- 
ning a s effects and his comments on the manuscript. 
SH thanks the HEPTOOLS Marie Curie RTN (contract 
number MRTN-CT-2006-035505) for an Early Stage Re- 
searcher position. TT thanks STFC (formerly PPARC) 
for an Advanced Fellowship. Financial support by MCNet 
(contract number MRTN-CT-2006-035606) and BMBF 
are acknowledged. 



A Alternative algorithm for phase 
space integration 

We explain in this section a method to integrate over the 
n-particle phase space, which was employed to cross-check 
the algorithm presented in Sec. We use an iterative 
approach to generate the event topology for the process 
PaPt — * Pi ■ ■ ■ Pn- For each step in the iteration we consider 
a2-> 2-scattering. Previous steps are taken into account 
by combining the particle momenta p a ,Pi ■ ■ - Pi into p ai 
and thereby considering the 2 — > 2-process p ai Pb PiPn- 



When denoting by 



and Si± the squared mass and 



squared transverse mass of the particle i, in the centre of 
mass frame of p ai b we obtain the integration boundaries 



1 



2 m„ 



(s ai b + Si - s n ) 



4 S aib 



where A 2 (s, si, S2) — (s — s\ — S2) 2 — 4siS2. 
sponding rapidity interval is fixed by 



The corre- 



V. 



„„„ = 1 ^ 1 + y/1 - S i± /ET^ 

2 l - ^l-^/f^J 



r 1 ± selec- 



and may be computed once k 2 ± is selected. The k 2 j 
tion is performed employing a divergence- free distribution, 
such as (k 2 jJ Q , where a > — 1. Since the above bound- 
aries are unambiguously determined, the n-particle phase 
space may be completely filled. 
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